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Abstract 

Increasing dispersal under range expansion increases invasion speed, 
which implies that a species needs to adapt more rapidly to newly 
experienced local conditions. However, due to iterated founder effects, local 
genetic diversity under range expansion is low. Evolvability (the evolution 
of mutation rates) has been reported to possibly be an adaptive trait itself. 
Thus, we expect that increased dispersal during range expansion may 
raise the evolvability of local adaptation, and thus increase the survival of 
expanding populations. We have studied this phenomenon with a spatially 
explicit individual-based metapopulation model of a sexually reproducing 
species with discrete generations, expanding into an elevational gradient. 
Our results show that evolvability is likely to evolve as a result of spatial 
variation experienced under range expansion. In addition, we show that 
different spatial phenomena associated with range expansion, in this case 
spatial sorting / kin selection and priority effects, can enforce each other. 
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Introduction 

Many species are currently expanding their ranges, both polewards and uphill, as a 
response to increasing global temperatures under climate change (Parmesan and Yohe 
2003; Chen et al. 2011). Range expansions are known to have profound effects on 
the genetic composition of populations, regarding both neutral and adaptive genetic 
diversity (Hewitt 1996; Thomas et al. 2001; Travis and Dytham 2002; Edmonds et 
al. 2004; Klopfstein et al. 2006; Phillips et al. 2006; Travis et al. 2007; Excoffier 
et al. 2009; Cobben et al. 2011; Cobben et al. 2012b). Traits that are to increase 
species dispersal capabilities and population growth rates are selected for under 
range expansions (Parmesan 2006; Burton et al. 2010; Phillips et al. 2010b; Hill et 
al. 2011; Shine et al. 2011). This may lead to the evolutionary increase of dispersal 
rate (Thomas et al. 2001; Travis and Dytham 2002; Kubisch et al. 2010; Henry et 
al. 2014), dispersal distance (Phillips et al. 2006) and effective fertility (Moreau et 
al. 2011) during periods of range expansion. In contrast, traits that affect the local 
adaptation of individuals may be affected by gene surfing (Cobben et al. 2012b), which 
is a consequence of the demographic processes occurring during colonization (Edmonds 
et al. 2004; Klopfstein et al. 2006), but the selective pressure on such traits is resulting 
only from the location of the individual and not the expansion itself. An increasing 
dispersal rate under range expansion will increase the invasion speed (Travis and 
Dytham 2002; Phillips et al. 2006) and has the implication that a species needs to be 
able to adapt to newly experienced local conditions more rapidly than before. However, 
Kubisch et al. (2013a) showed that high dispersal rates prevent the evolution of local 
adaptation. Now there is both theoretical and empirical evidence that evolvability 
can be adaptive itself under conditions that require an increased rate of adaptation, 
e.g. under increasing environmental stochasticity and stress (Leigh Jr 1970; Ishii et al. 
1989; Sniegowski et al. 1997; Kashtan et al. 2007; Lee and Gelembiuk 2008). 

Evolvability generally refers to the ability of populations to evolve in an adaptive way 
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(Brookfield 2009) and is associated with varying environments (Kashtan et al. 2007). 
Many studies have speculated and shown the increase of evolvability under temporally 
varying environmental conditions (Leigh Jr 1970; Ishii et al. 1989; Sniegowski et al. 
1997; Kashtan et al. 2007). Under such conditions, a high mutation rate serves to 
replenish the local genetic variation, which has beforehand been strongly reduced 
under differing selection regimes (Brookfield 2009). With a high frequency of local 
disturbances, the mutation rate is expected to be maintained at a constantly high 
level (Ishii et al. 1989; Earl and Deem 2004). In contrast, when we take a spatial 
perspective, spatial heterogeneity may lead to differing selection pressures between 
habitat patches. High levels of dispersal can thus result in individuals experiencing 
a varying environment in space rather than time (Lee and Gelembiuk 2008). In a 
scenario of range expansion, increasing relatedness between individuals reduces local 
genetic diversity (Kubisch et al. 2013b), while dispersal rates increase. Under such 
low local genetic diversity, this high dispersal rate should then be accompanied by high 
mutation rates, i.e. high evolvability, to increase local adaptation and survival. 

In this study, we have sought to elaborate on these thoughts and used a simulation 
model to investigate the interplay between the evolution of dispersal rate and the 
evolution of evolvability under range expansion across a spatial gradient in altitude. 



The Model 

We are using a spatially explicit individual-based metapopulation model of a sexually 
reproducing species with discrete generations distributed along an elevational gradient. 
The basic model has already been successfully applied in theoretical studies, mainly 
focused on dispersal evolution (Travis et al. 1999; Kubisch et al. 2010; Fronhofer et 
al. 2011; Kubisch et al. 2013a; Kubisch et al. 2013b) and was parameterized using 
empirical data (Poethke et al. 1996; Amler et al. 1999). For the current study we allow 
evolvability of local adaptation to be adaptive, and investigate the interplay with the 
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evolution of dispersal rate. 

Landscape 

The simulated landscape consists of 250 columns (x-dimension) of 20 patches each 
(y-dimension) . We assume wrapped borders, building a torus. Hence, if an individual 
leaves the world in y-direction during dispersal, it will reenter the simulated world on 
the opposite side. However, if it leaves the world in the x-direction, it is lost from the 
simulation. While most studies investigating range expansions across an environmental 
gradient, for the sake of simplicity focus on a single parameter changing along space, 
a typical elevational gradient is known to involve both a decreasing temperature and 
increasing habitat fragmentation (Korner and Paulsen 2004). Thus, in our model 
firstly every column of patches (x-position) is characterized by its specific abiotic 
habitat conditions r x . Throughout this manuscript, r x will be interpreted as 'mean 
temperature'. This mean local temperature is used for the determination of local 
adaptation of individuals. To simulate a large-scale habitat gradient, x changes linearly 
from t x= \ = 0 to ^=250 = 0 along the x-dimension, i.e. by A TX = 0.04 when 
moving one step in x-direction. Secondly, to account for habitat fragmentation in terms 
of patch isolation, each x-position is characterized by a certain degree of dispersal 
mortality. The probability to die upon emigration \x changes linearly from \x x= \ = 0 
to fi x =25o = 1 along the x-dimension. 

Population dynamics and survival of offspring 

Local populations are composed of individuals that are characterized by several 
traits: 1) their sex, 2) two alleles at locus 1 coding for the individuals emigration 
propensity, 3) two alleles at locus 2 coding for the individual's 'habitat preference', 
i.e. the environmental conditions (temperature r) under which the individual survives 
best (see below for details), and 4) another 2 alleles at locus 3 coding for the mutation 
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probability of the alleles at locus 2, i.e. the evolvability of habitat preference (see below 
under Genetics). 

Local population dynamics follow the time-discrete BevertonHolt model (Beverton 
and Holt 1957). Each individual female in patch x,y is therefore assigned a random 
male from the same habitat patch (males can potentially mate several times) and gives 
birth to a number of offspring drawn from a Poisson distribution with mean A XjVt t. 
The offspring's sex is chosen at random. Density-dependent survival probability si of 
offspring due to competition is calculated as: 



Finally, the surviving offspring experience a further density-independent mortality 
risk (1 — S2) that depends on the matching of their genetically determined optimal 
temperature (r opt ) to the temperature conditions in patch x, y (r x ) according to the 
following equation: 



Density-independent survival depends on adaptation to local conditions, i.e. on the 
difference between the genetically encoded optimal temperature r opt of an individual 
and local temperature conditions t x . r\ describes the niche width or 'tolerance' of 
the species. We performed simulations for the species with a niche width of rj = 
0.5, equivalent to a decrease of survival probability of about 0.02 when dispersing 
one patch away from the optimal habitat. By using this approach we assume that 
density-dependent mortality (1 — si) acts before mortality due to maladaptation to 
local conditions (1 — s 2 ). 

Individual surviving offspring disperse with probability d that is determined by their 
dispersal locus (see below). If an individual emigrates it dies with probability //, which 



1 



(1) 



si = 



l + ^-Nx,y,t 
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is calculated as the arithmetic mean of the dispersal mortality values of its natal and 
its target patch, respectively. This mortality accounts for various costs that may be 
associated with dispersal in real populations, like fertility reduction or predation risk 
(Bonte et al. 2012). We assume nearest-neighbor dispersal, i.e. successful dispersers 
settle randomly in one of the eight surrounding habitat patches. 

Genetics 

As mentioned above, every individual carries three unlinked, diploid loci coding for 
its emigration probability, its habitat preference (optimum temperature), and the 
mutation rate of the optimum temperature alleles (evolvability of habitat preference), 
respectively. The phenotype of an individual is determined by calculating the 
arithmetic means of the two corresponding alleles. Hence, dispersal probability d is 
given by d = ( wit h i dl an d l d2 giving the 2 'values' of the two dispersal 

alleles), optimal temperature T opt is calculated as r opt = la ' 1 + la ' 2 (with l a> i and l a ^ 
giving the 'values' of the two adaptation alleles), and similarly the mutation rate of 
optimal temperature m{r opt ) = 10~ ex ' p (with exp = le ' 1 ~^ le ' 2 ; and and / 6j2 the 
'values' of the two evolvability alleles). At each of the three loci, newborn individuals 
inherit alleles, randomly chosen, from the corresponding loci of each of their parents. 
During transition from one generation to the next an allele may mutate. Alleles at 
the dispersal locus and the evolvability locus mutate with a probability of m = 
10~ 4 . Alleles at the adaptation locus however, mutate with the probability m(r opt ) 
given by the value based on its two alleles at the evolvability locus as elaborated 
above. Mutations are simulated by adding a random number drawn from a Gaussian 
distribution with mean 0 and standard deviation 0.2 (in case of the alleles coding for 
optimal temperature the standard deviation is 0.5) to the value calculated from the 
mean of the inherited alleles. 
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Simulation Experiments 

Simulations were initialized with a 'native area' (from x = 1 to x — 50) from where 
the species was able to colonize the world, while the rest of the world was initially kept 
free of individuals. Upon initialization, dispersal alleles (ld,i) were randomly drawn 
from the interval 0 < l^i < 1, and evolvability alleles were set to 4, added a 
Gaussian random number with mean zero and standard deviation one. Populations 
were initialized with K locally optimally adapted individuals, i.e. preference alleles 
were initialized according to the local temperature r x . However, to account for some 
standing genetic variation we also added to every respective optimal temperature allele 
a Gaussian random number with mean zero and standard deviation 0.5. We performed 
100 replicate simulations, which all covered a time span of 60, 000 generations. To 
establish equilibrium conditions, individuals were confined to their native area during 
the first 1, 000 generations. After this burn-in period, the species was allowed to pass 
the x = 50 border. Table 1 summarizes all relevant model parameters, their meanings 
and the standard values used for the simulations. 

Analysis 

The individual phenotypes for the three traits were documented in time and space 
throughout the simulations. Genetic diversity was calculated as the variance in 
allelic values at the adaptation locus per x-position. The marginal values of dispersal 
propensity and evolvability were calculated as the arithmetic mean of all individual 
values at the range border, i.e., the mean of all patches in the y dimension of the last 
five x-positions of the gradient counted from the most-forward occupied patch (which is 
time dependent). 
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Table 1: Used parameter values. 



parameter/variable value meaning 

individual variables: 

ld,i,ld,2 evolving alleles coding for the dispersal propensity 





evolving 


alleles coding for the optimal temperature 




evolving 


alleles coding for the mutation rate of the 






optimal temperature 


simulation parameters: 






K 


100 


carrying capacity 


A 


2 


per capita growth rate 


e 


0.05 


local extinction probability 


m 


10" 4 


mutation rate for dispersal and evolvability 






alleles 




[0..1] 


local dispersal mortality 


r x 


[0..10] 


local temperature 


V 


0.5 


niche width 


■Kmax 


250 


extent of simulated landscape in x-direction 


Umax 




extent of simulated landscape in y-direction 



Results 

After the burn-in phase, the dispersal propensity d in the core area, so under low 
dispersal mortality, was on average approximately 0.35. Maximum population 
density was around 0.85 here, accompanied by a high level of local adaptation S2 (i.e. 
adaptation-dependent offspring survival probability) close to one (not shown), the 
mutation rate of the optimum temperature m(r opt ) remained low, between 10~ 4 and 
10" 5 . 

Under range expansion, the dispersal rate d increased in the populations at the 
expansion front, and in time decreased again as these populations got older and 
adapted locally (Figure IB). During range expansion the maximum established 
dispersal rate d was approximately 0.25 (Figure IB). As the species expanded its 
range further across the gradient of increasing dispersal mortality (and increasing 
temperature), the maximum established dispersal rates at the range front decreased 
(Figure IB). The range border established at an average dispersal mortality of 0.8, and 



8 



Downloaded from http://biorxiv.org/ on September 18, 2014 



Cobben & Kubisch 




~i 1 1 1 1 n 1 1 1 1 — 

0 50 100 150 200 0 50 100 150 200 

spatial location 



Figure 1: The average values over 100 simulations during and after range expansion 
across the gradient (horizontal axis) in time (gray scaling from light to dark, as time 
proceeds, which is given in a sequence of generations 1000, 1250, 1500, 2000, 5000, 
10000, 60000) of A. population density, B. emigration propensity, C. evolvability, 
i.e. mutation rate of the adaptation alleles, and D. genetic diversity measured as the 
variance in adaptation alleles. Shown are mean values. For reasons of clarity, a moving 
average with a window size of 20 has been applied (data were present in 10-generation 
intervals). 

a temperature value of 8. After reaching spatial equilibrium, the spatial distribution 
of dispersal rates showed a decaying exponential trend. The mutation rate of the 
temperature optimum m(r op t) showed no spatial pattern before range expansion, 
but rapidly increased once the populations invaded into the landscape. After 5, 000 
generations of expansion, this mutation rate was on average almost seven times 
higher at the range margin than in the core (Figure 1C). As a result of this increased 
mutation rate, genetic diversity was also increased at the range front, compared to 
regions, which have been populated for longer time (Figure ID). Diversity in the initial 
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d 




50 200 1000 5000 

generations 



Figure 2: The average values over 100 simulations at the range border of A. the 
emigration propensity and B. the mutation rate in time (note the logarithmic scale). 
Shown are mean values (solid black lines) and 25%- and 75%-quartiles (dashed grey 
lines). For reasons of clarity, a moving average with a window size of 20 has been 
applied. For details see main text. 

core area, however, stayed at fairly high levels. Similar to mean dispersal propensity, 
mutation rates decreased again when the populations were getting older (Figure 1C). 
However, the time lag between the local decrease of the dispersal rate and the decrease 
of the mutation rate m(r opt ) was steadily increasing across space. At the range border 
it took 45, 000 generations for the mutation rate to decrease to equilibrium values 
after the dispersal rate had decreased (Figure 2A,B). Note also that variability in 
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mutation rates was very high. While the maximum average mutation rate was found 
to be around 2 • 1CT 4 , the third quartile was measured at about 1.3 • 10~ 3 . The local 
level of adaptation S2 was close to one throughout the simulation time and across the 
complete species range. 

Discussion 

In this study we investigated, whether an increase in dispersal rate under range 
expansion can lead to the evolution of higher evolvability, enabling species to more 
rapidly adapt to local conditions with an increasing invasion speed. We found an 
increase of evolvability, which is dependent on the local level of genetic diversity, with 
high mutation rates evolving under high rates of dispersal. Evolved high mutation 
rates took an extensive period of time to return to lower mutation rates after the local 
dispersal rates decreased again, due to a priority effect (De Meester et al. 2002; Urban 
and De Meester 2009). 

This modeling study shows for the first time that evolvability can evolve as a result of 
spatial variation, experienced under range expansion. This high evolvability increases 
genetic diversity and thus adaptive potential in newly colonized areas. In addition, 
we show that different spatial phenomena associated with range expansion, in this 
case spatial sorting / kin competition and priority effects, can enforce each other. 
The genetic signature of spatial sorting is then extended in longevity due to a priority 
effect, where high local adaptation further delays the establishment of lower dispersal 
rates. This may have implications for the interpretation of field data. 

During range expansion the dispersal rate was showing a clear signal of spatial sorting 
(Shine et al. 2011) and kin competition (Kubisch et al. 2013b), with good dispersers 
gathering at the expanding wave front (Phillips et al. 2010a). The immigration of 
different individuals is expected to maintain a high local level of genetic variation 
(Holt and Barfield 2011), from which one would expect high levels of dispersal to 
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be accompanied by a low local mutation rate and we indeed find this pattern in the 
range core throughout the simulation. At the margin, however, relatedness amongst 
individuals increases at an advancing range front (Kubisch et al. 2013b), reducing 
both local genetic diversity and the diversity of immigrants. Under these conditions 
an increase of the evolvability of local adaptation evolved, which compensated for the 
experienced spatial variation in local temperature. Both high rate signals of dispersal 
and mutation rate disappeared with time. At the range border, individuals with a high 
dispersal rate were locally well-adapted as a result of the high mutation rate. While 
a lower dispersal rate is beneficial under conditions of high dispersal mortality, slow 
dispersers took a long time to reach the area (genetic signature of range expansion, 
Phillips et al. 2010b). They were in addition hindered by the high dispersal mortality, 
and on top of that needed to compete with better locally adapted individuals (priority 
effect, De Meester et al. 2002). This has implications for field work, where the cause 
of an observed local high dispersal rate requires careful interpretation, as it can be 
the result of natural selection, spatial selection, spatial disequilibrium, or priority 
effects. The dispersal rate decreased first and was only after an extensive time lag 
of ten thousands of generations followed by a decrease in mutation rate. This again 
was caused by a priority effect (De Meester et al. 2002), where a high level of local 
adaptation prevents the establishment of individuals with a lower mutation rate. 
We speculate that this second time lag was larger than the first because selection on 
evolvability is weaker than on dispersal, with dispersal having more direct and drastic 
effects on population dynamics. However, this remains to be investigated. 

Holt and Barfield (2011) investigated niche evolution at species range margins and 
found that local evolution was hampered when source populations of immigrating 
individuals were at low density, as a result of the stochastic processes in such 
populations (Pearson et al. 2009; Bridle et al. 2010; Turner and Wong 2010). The 
likelihood of observing niche evolution was further affected by the mutation rate, where 
dispersal limited local evolution in the sink population under a higher mutation rate, 
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because of the increased numbers of maladapted individuals from the source (Holt 
and Barfield 2011). They did, however, not allow the joint evolution of mutation rate 
and dispersal rate, but instead used fixed rates. Our results can be affected by the 
used genetic architecture, where linkage between traits (Blows and Hoffmann 2005; 
Hellmann and Pineda-Krch 2007), polygeny, and the magnitude of mutations can be of 
importance in range dynamics (Kawecki 2000; Kawecki 2008; Walsh and Blows 2009; 
Gomulkiewicz et al. 2010; Kimbrell 2010). Increased evolvability can be modelled 
in different ways, e.g. an increased magnitude of the phenotypic effect of mutations 
(Griswold 2006), the evolution of modularity (Kashtan et al. 2009), or the evolution of 
generalism or plasticity (Lee and Gelembiuk 2008; Chevin and Lande 2011), where we 
have restricted evolvability to the mutation rate. Although empirical studies support 
the decision to focus on mutation rate (Sniegowski et al. 2000; Earl and Deem 2004), 
this can have an effect on our results. 

There is an ever-expanding pool of literature discussing the ecological and evolutionary 
dynamics of dispersal in the formation of species ranges (reviewed in Kubisch et al. 
2014). The evolution of dispersal has been shown to increase invasion speeds (Thomas 
et al. 2001; Travis and Dytham 2002; Phillips et al. 2010a), affect the fate of neutral 
mutations (Travis et al. 2010), as well as the level of local adaptation (Kubisch et al. 
2013a; Bourne et al. 2014), and local population dynamics (Travis et al. 2007; Burton 
et al. 2010), and in addition cause strong patterns of spatial disequilibrium (Ibrahim 
et al. 1996; Phillips et al. 2010b). While individual-based models have recently largely 
extended our theoretical knowledge of interactions and evolution of traits during range 
expansion, empirical data have been restricted to a few well-known cases (Thomas et 
al. 2001; Phillips et al. 2006; Moreau et al. 2011). Increasing ecological realism in our 
models (Cobben et al. 2011; Cobben et al. 2012a; Cobben et al. 2012b; Bocedi et al. 
2014) might improve the predictability of theoretical phenomena and support field 
studies. These are, however, always constrained by the required temporal and spatial 
scales, which are particularly restrictive in terrestrial systems. 
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